home *** CD-ROM | disk | FTP | other *** search
/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / gaussian.t < prev    next >
Encoding:
Text File  |  1997-04-18  |  2.8 KB  |  163 lines
  1. %
  2. % "gaussian.t" performs gaussian elimination on 
  3. % an N x N matrix and solves for x in:
  4. %
  5. %                   A x = b
  6. %
  7. %   Sample program for the T Interpreter by:
  8. %
  9. %   Stephen R. Schmitt
  10. %   962 Depot Road
  11. %   Boxborough, MA 01719
  12. %
  13.  
  14. const N : int := 3
  15.  
  16. var a : array[N,N+1] of real
  17. var x : array[N] of real
  18.  
  19. program
  20.  
  21.     var i, j : int
  22.  
  23.     a[0,0] :=  1        % a00
  24.     a[0,1] :=  3        % a01
  25.     a[0,2] := -4        % a02 
  26.     a[0,3] :=  8        % b0  
  27.  
  28.     a[1,0] :=  1        % a10
  29.     a[1,1] :=  1        % a11
  30.     a[1,2] := -2        % a12
  31.     a[1,3] :=  2        % b1 
  32.  
  33.     a[2,0] := -1        % a20
  34.     a[2,1] := -2        % a21
  35.     a[2,2] :=  5        % a22
  36.     a[2,3] := -1        % b2 
  37.  
  38.     x[0] := 0           % initialize x 
  39.     x[1] := 0
  40.     x[2] := 0
  41.  
  42.     if eliminate then   % lower diagonal elements
  43.         
  44.         substitute      % and solve for x
  45.  
  46.         put "solution:"
  47.         put ""
  48.         for i := 0 ... N-1 do
  49.             put "    x[", i, "] = ", x[i] 
  50.         end for
  51.         
  52.     else
  53.     
  54.         put "singular matrix!"
  55.         
  56.     end if
  57.  
  58. end program
  59.  
  60. %
  61. %   "eliminate" converts A to upper diagonal form
  62. %
  63. %   returns:  nothing
  64. %
  65. function eliminate : boolean
  66.  
  67.     var i, j, k, max_row : int
  68.     var t : real
  69.     var ok : boolean := true
  70.  
  71.     for i := 0 ... N-1 do
  72.  
  73.         %
  74.         % find row with maximum in column i
  75.         %
  76.         max_row := i
  77.  
  78.         for j := i+1 ... N-1 do
  79.  
  80.             if abs( a[j,i] ) > abs( a[max_row,i] ) then
  81.  
  82.                 max_row := j
  83.  
  84.             end if
  85.  
  86.         end for
  87.  
  88.         %
  89.         % swap max row with row i of A and b
  90.         %
  91.         for k := i ... N do
  92.   
  93.             t            := a[i,k]
  94.             a[i,k]       := a[max_row,k]
  95.             a[max_row,k] := t
  96.       
  97.         end for
  98.  
  99.         %
  100.         % eliminate lower diagonal elements
  101.         %
  102.         for j := i+1 ... N-1 do
  103.  
  104.             for decreasing k := N ... i do
  105.  
  106.                 if a[i,i] = 0.0 then
  107.                     ok := false
  108.                 else
  109.                     a[j,k] := a[j,k] - a[i,k] * a[j,i] / a[i,i]
  110.                 end if
  111.  
  112.             end for
  113.  
  114.         end for
  115.  
  116.     end for
  117.  
  118.     return ok
  119.  
  120. end function
  121.  
  122. %
  123. %   "substitute" computes the values of x starting from the bottom
  124. %
  125. %   returns:  nothing
  126. %
  127. procedure substitute
  128.  
  129.     var j, k : int
  130.     var t : real
  131.  
  132.     for decreasing j := N-1 ... 0 do
  133.  
  134.         t := 0.0
  135.  
  136.         for k := j+1 ... N-1 do
  137.  
  138.             t := t + a[j,k] * x[k]
  139.  
  140.         end for
  141.  
  142.         x[j] := ( a[j,N] - t ) / a[j,j]
  143.  
  144.     end for
  145.  
  146. end procedure
  147.  
  148. %
  149. %   "abs" returns the absolute value of a real argument
  150. %
  151. %   returns:  real number
  152. %
  153. function abs( x : real ) : real
  154.  
  155.     if x < 0 then
  156.  
  157.         x := -x
  158.  
  159.     end if
  160.  
  161.     return x
  162.  
  163. end function